How self-organized criticality works: A unified mean-field picture 
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We present a unified mean-field theory, based on the single site approximation to the master- 
equation, for stochastic self-organized critical models. In particular, we analyze in detail the proper- 
ties of sandpile and forest-fire (FF) models. In analogy with other non-equilibrium critical phenom- 
ena, we identify the order parameter with the density of "active" sites and the control parameters 
with the driving rates. Depending on the values of the control parameters, the system is shown to 
reach a subcritical (absorbing) or super-critical (active) stationary state. Criticality is analyzed in 
terms of the singularities of the zero-field susceptibility. In the limit of vanishing control parameters, 
the stationary state displays scaling characteristic of self-organized criticality (SOC). We show that 
this limit corresponds to the breakdown of space-time locality in the dynamical rules of the models. 
We define a complete set of critical exponents, describing the scaling of order parameter, response 
functions, susceptibility and correlation length in the subcritical and supercritical states. In the 
subcritical state, the response of the system to small perturbations takes place in avalanches. We 
analyze their scaling behavior in relation with branching processes. In sandpile models because of 
conservation laws, a critical exponents subset displays mean-field values {v = 1/2 and 7 = 1) in any 
dimensions. We treat bulk and boundary dissipation and introduce a new critical exponent relating 
dissipation and finite size efi'ects. We present numerical simulations that confirm our results. In the 
case of the forest-fire model, our approach can distinguish between different regimes (SOC-FF and 
deterministic FF) studied in the literature and determine the full spectrum of critical exponents. 

PACS numbers; 64.60.Lx, 05.40. -f-j, 05.70.Ln 



I. INTRODUCTION 

After ten years of research and countless papers the 
precise significance of self-organized criticality (SOC) |^] 
is still controversial. Originally, SOC was presented as a 
general theory to understand fractals and 1// noise as the 
natural outcome of the dynamical evolution of systems 
with many coupled degrees of freedom. Irreversible dy- 
namics would generate a self-organization of the system 
into a critical state, without the fine tuning of external 
parameters. The SOC idea was illustrated by computer 
models in which a slow external driving leads to a sta- 
tionary state with avalanches of widely distributed ampli- 
tude [|l|. This proposal stimulated a cascade of research 
activity in experiments, theory and simulations. While 
the explanation presented in Ref. [Q about the origin of 
scaling in nature appears now too simplistic, SOC gave 
a formidable input to the study of slowly driven systems 
and avalanche phenomena. 

Avalanche behavior was experimentally observed in a 
variety of phenomena ranging from magnetic systems 
(the Barkhausen effect) Q and flux lines in high-Tc su- 
perconductors 1^], fluid flow through porous media [Q, 
microfracturing processes earthquakes and lung 
inflation M. In addition, SOC ideas stimulated a great 



interest in granular matter [^j , although it was soon real- 
ized that the concept was hardly applicable there, apart 
from the academic example of a ricepile |^ . All the men- 
tioned experiments share with SOC models the slow ex- 
ternal driving and the avalanche response, but it is un- 
clear whether self-organization as described in Ref. 0] 
plays any role in there. To answer this question it would 
be necessary to better understand what determines the 
appearance of scaling in SOC models and driven systems 
in general. 

The idea of a critical point without flne tuning of ex- 
ternal parameters is very appealing because it opposes 
the standard picture of equilibrium critical phenomena. 
The concept of "spontaneous" criticality, as it has been 
discussed in the literature, presents, however, several am- 
biguities. It has been pointed out by several authors 
that the driving rate is a parameter that has to be fine 
tuned to zero in order to observe criticality JlO|-p^ . This 
fact poses no problems to computer simulations, where 
an infinite timescale separation can easily be enforced, 
but it is crucial in experiments where the driving rate 
is always non-zero. The second ambiguity is mostly a 
language problem: calling "self-organization" the evolu- 
tion towards a stationary state can be misleading. Any 
non-equilibrium system poised at its "fine tuned" critical 
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point, when started from a generic configuration, evolves 
towards the critical stationary state, thus building up 
correlations and scaling. We would not describe this pro- 
cess as self-organization. These ambiguities in the defi- 
nition of SOC have hindered the formulation of precise 
relations with other non-equilibrium critical phenomena 

0- 

In the past years, several attempts have been made to 
find a general mechanism to describe SOC models. Sor- 
nette proposed that SOC was due to a non-linear feed- 
back between order parameter and control parameter, 
leading to the self-tuning of the latter to the critical point 
. Later it was pointed out that this mechanism could 
lead also to a first-order transition rather than a criti- 
cal point Ip^ . In a recent paper, it was claimed instead 
that SOC corresponds to the tuning to zero of the order 
parameter of an ordinary critical phenomenon . Our 
analysis shows that the situation is simpler: criticality 
arises from the fine tuning to zero of one or more con- 
trol parameters (driving rate, dissipation) and there is 
no coupling between control and order parameters [ p^ . 
The incorrect identification of control and order param- 
eters is at the basis of many misconceptions about SOC 
phenomena, as we will discuss in the following. 

Many theoretical methods have been used in the anal- 
ysis of SOC models. Few rigorous milestones can be 
found in the activity of Dhar and coauthors and 
in the Dubna group works [|9| . Flory and Langevin- 
type approaches pl|-p3[ have been used on a phenomeno- 
logical basis. More recently, real space renormalization 
group provided good estimates of the avalanche expo- 
nents [|l2|,^. Despite their richness, however, all these 
approaches are focused on the critical avalanche behav- 
ior and the external driving does not play any role; i.e. 
the system is studied in the infinite timescale separation 
regime. Furthermore, many of these attempts are con- 
ceived ad hoc for particular models and do not provide a 
general conceptual framework for SOC phenomena. 

The first step towards a comprehensive theoretical un- 
derstanding of SOC is provided by mean-field (MF) the- 
ory which gives insight into the fundamental physical 
mechanisms of the problem and a reference language. 
It provides a feasible treatment to nonequilibrium and 
complex problems (often the only one), and can be used 
as a starting point for more sophisticated calculations. 
Whereas many numerical and analytical approaches get 
harder as the dimensionality increases, MF theory im- 
proves and despite the crude approximations it usually 
gives correct qualitative predictions for the phase dia- 
grams of high-dimensional systems. Finally, MF theory 
highlights the importance of symmetries and conserva- 
tion laws. 

A vast activity concerning MF theory of SOC mod- 
els can be found in the literature. Exponents describing 
avalanche distributions and propagation have been com- 
puted in several ways: solving infinite-range p^ ], Bethe 
lattice and random neighbor p6|-p8[| models and by 
mapping the dynamics into a branching process [p9|-^ . 



Self-consistent MF approximations for the sand height 
distribution have also been used [ p0|j34| . Other MF ap- 
proaches use analogies with equilibrium critical phenom- 
ena |35y36|] , leading sometimes to incorrect predictions as 
we will discuss in the following. In summary, the MF ap- 
proach to SOC systems is composed of a number of stud- 
ies of specific models but a comprehensive understanding 
of the phenomenon is missing. 

Here, we present a unified MF description of SOC mod- 
els using the formalism developed for non-equilibrium 
critical phenomena with steady states. We use a single 
site approximation to the master-equation and we en- 
force conservation laws by effective parameters and con- 
straints. We concentrate on models driven by stochastic 
noise, such as the sandpile 0] and the forest-fire (FF) 
prf . In order to write a master equation, we consider fi- 
nite values of the driving rates, since only in this case the 
dynamical rules are local in space and time. Our analy- 
sis shows that criticality in these models corresponds to 
the limit in which the dynamical rules become non-local. 
Non-locality is implicitly enforced in computer simula- 
tions, where the evolution of a single site depends on the 
state of the entire system. This fact is particularly evi- 
dent in extremal models in quenched disorder, where the 
dynamics proceeds by a global minimum search. Also 
in that case, to write a local master-equation one has to 
introduce a non-zero driving rate, but the driving mech- 
anism differs from the one we discuss here and will be 
described in a forthcoming paper ]3S| ] . 

The present approach allows us to identify control and 
order parameters of SOC models and to clarify the re- 
lations with other non-equilibrium critical phenomena 
p3[ . In particular, we show that SOC models have close 
similarities with non-equilibrium cellular automata with 
many absorbing states |pO|-^. The major difference is 
that in SOC models, the control parameters have to be 
tuned to zero to reach criticality. As we discussed before, 
this limit corresponds to the breakdown of locality in the 
dynamical rules of the model and hence to the onset of 
long range correlation in the dynamical response. While 
apparently there is no mathematical difference between 
tuning a parameter to zero or to a non-zero value, the 
physical differences are quite important 0]. In the first 
case, changing the value of the control parameter by a 
given factor still keeps the system close to the critical 
point. This is not the case in ordinary phase transitions, 
where, doubling the value of the temperature, the sys- 
tem looses completely the critical properties. Moreover, 
in order for SOC model to be defined, the control pa- 
rameter (i.e. the driving rate) should always be non-zero 
and the critical point can only be reached through a limit 
process. 

The present MF theory can be applied to any stochas- 
tic cellular automaton and therefore provides a unified 
description for the ensemble of SOC models and other 
related non-equilibrium critical systems, such as contact 
processes and cellular automata with absorbing states 
pO|-^. Moreover, it serves to emphasize the differences 
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between different models and between different regimes 
in the same model. Our analysis points also out the in- 
consistencies contained in earlier MF approaches [ ^J36| 
which led to a misleading characterization of the model. 
We identify the sub-critical and super-critical states of 
SOC models and discuss the different ways in which crit- 
icality can be reached. We describe the avalanche behav- 
ior characteristic of these models in terms of response 
functions and we study the effect of perturbations on the 
stationary state. We introduce a full set of critical expo- 
nents, describing the response at the critical point and 
the scaling close to the critical point in the sub-critical 
and super-critical states. In the case of sandpile models 
a subset of exponents is found to have mean-field values 
in any dimension. The reason for this behavior, which 
is confirmed by numerical simulations, is ascribed to the 
presence of conservation laws in the dynamics. 

The paper is organized as follows: In Section II we in- 
troduce the models. In Section III we review the dynamic 
mean-field approximation to the master equation. Sec- 
tion IV contains the mean-field theory for the sandpile 
model and discusses some issues related to conservation. 
In Section V we report the mean-field analysis for the FF 
model and Section VI is devoted to a general discussion. 
A brief report of these results has appeared in Ref. |13|. 



II. THE MODELS 



A rapid look at the SOC literature discourages every 
newcomer in the field. In less than ten years more than 
two thousand papers have been published, and compre- 
hensive reviews are not yet appeared (a valuable effort 
in this direction can be found in Ref This is due 

to the lack of a general understanding which would pro- 
vide the framework to order the huge amount of infor- 
mations about SOC. In particular, we were spectators of 
an hectic activity in numerical simulations, with the in- 
troduction of a multitude of different models. A closer 
look at the literature reveals that the number of original 
models can be greatly reduced by noting that most of 
them are variations of prototype models. Using a more 
"Draconian" approach, we can distinguish just two main 
families of SOC models. The first is represented by the 
stochastic SOC models such as the sandpile or forest-fire 
model, in which the self-organization process is the out- 
put of a stochastic dynamics. The second family groups 
together the so-called "extremal" or "quenched" models 
p8[ , which are defined by a deterministic dynamics in a 
random environment. Examples of the latter family are 
the Invasion Percolation (IP) and the Bak-Sneppen 
(BS) iQ models. In this paper we discuss stochastic 
models, but work is in progress to extend the present 
analysis to systems driven by an extremal dynamics. 



A. Sandpile models 

Sandpile models are cellular automata (CA) with an 
integer (or continuous) variable zi (energy) defined in a 
d— dimensional lattice. At each time step an energy grain 
is added to a randomly chosen site, until the energy of a 
site reaches a threshold Zc- When this happens the site 
relaxes 

Z'l > Z{ Zq (1) 

and the energy is transferred to the nearest neighbors 
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(2) 



The relaxation of a site can induce nearest neighbor sites 
to relax on their turn, i.e. they exceed the threshold 
because of the energy received. New active sites can 
generate other relaxations and so on, eventually giving 
rise to an avalanche. For conservative models the trans- 
ferred energy equals the energy lost by the relaxing site 
i^Uj = Zc), at least on average. Usually, the only form 
of dissipation occurs at the boundary, from which en- 
ergy can leave the system. It is worth to remark that 
during the avalanches, the energy input stops until the 
system is again in equilibrium and no active sites are 
presents. This corresponds to an infinite timescale sepa- 
ration. With these conditions the system reaches a sta- 
tionary state characterized by avalanches whose sizes s 
are distributed as a power law p],^^- 



P{s) 



(3) 



The model originally introduced by Bak, Tang and 
Wiesenfeld (BTW) [§ is a discrete automaton in which 
Zc — 2d and yj = 1. An interesting variation of the 
original sandpile is the three-states Manna model p9[ . 
In this automaton the critical threshold is Zc = 2 in- 
dependently on the dimensionality d and if a relaxation 
(toppling) takes place, the energy is distributed into two 
randomly chosen nearest neighbor sites. Other variations 
in which part of the energy is kept by the relaxing site 
can also be considered as well as directed models in which 
energy is transferred along a preferential direction p5[ . 
Finally, sandpile models that include a relaxation dynam- 
ics where part of the energy is dissipated have been con- 
sidered ll^. These models can be characterized by the 
fraction of energy that disappears from the system dur- 
ing each relaxation process. When a global dissipation 
is present (energy is lost on average), the critical behav- 
ior is destroyed and a characteristic length is introduced. 
These numerical evidences suggest that conservation is 
necessary to obtain criticality. 

As we discussed above, sandpile models are driven by 
adding a single energy grain on a randomly chosen site, 
when no active site is present. In this way, avalanches are 
instantaneous with respect to the driving timescale. This 
rule is very naturally implemented in a computer algo- 
rithm, which can handle at the same time the two differ- 
ent timescales. This, however, corresponds to a nonlocal 
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interaction in which the site dynamical evolution depends 
upon the whole system configuration. This non-local in- 
teraction is hard to describe, and in order to perform an 
analytical description we have to fix a reference timescale, 
as for example the single site relaxation step, and mea- 
sure the driving rate on that scale. For this reason, we 
consider a generalized sandpile model |^ , that includes 
a non vanishing driving rate by introducing the proba- 
bility h per unit time that a site will receive a grain of 
energy. Energy is distributed homogeneously and the to- 
tal energy flux is given by Ji„ = hL"^. The parameter 
h sets the driving timescale or equivalently the typical 
waiting time between different avalanches as ~ 1/h. 
In the limit h 0, we recover the slow driving limit; i.e. 
during an avalanche the system does not receive energy. 
This formulation of the dynamics has the advantage to be 
local in space and time. The state of a single site depends 
only on the state of the site itself and its nearest-neighbor 
sites at the previous time step, through a transition prob- 
ability that is given by the reaction and driving rates. It 
will be also convenient in the ensuing analysis to group 
the possible states of a site in three classes: active when 
z > Zc, critical when z = Zc — 1 and stable for z < Zc — 1. 



B. Forest-Fire model 

The first example of a stochastic SOC model without 
conservation is in the forest fire (FF) model |M. The 
model has been first introduced by Bak et al. pa as an 
example of SOC, and has been then modified by Drossel 
and Schwabl ||5^. The model is defined on a lattice in 
which each site can be empty, occupied by a green tree or 
by a burning tree. Burning trees turns to ashes with uni- 
tary rate and set fire to the neirest neighbor trees. The 
model was first studied in the case of a small tree growth 
rate p and in the absence of a spontaneous ignition of 
fires. In d = 2 the system reaches a dynamical state in 
which fire fronts propagate with trivial scaling proper- 
tics fs^ . Only recently new large scale simulations have 
shown that for d > 2 anomalous scaling laws occur . 
A more interesting situation appears when a very small 
rate for spontaneous fire ignition / (lightning probabil- 
ity) is introduced in the automaton dynamics. The sys- 
tem shows scaling behavior with a diverging characteris- 
tic length in the limit f /p and p ^ and the activity 
occurs in bursts of fire spreading (avalanches) whose dis- 
tribution follows a power law behavior P(s) ~ s~'^ . In 
the FF model the driving rates are explicitly defined by / 
and p, and the dynamical rules are thus local. However, 
in numerical simulations the two driving field are implic- 
itly set to zero by the condition that trees growth and fire 
ignition occur only when the system does not show active 
sites. Only the ratio 9 = f /p is quantitatively defined by 
the relative probability of tree growth with respect to fire 
ignition events. Also in this case numerical simulations 
are done in the infinite timescale separation limit, which 



corresponds to the subcritical state of the system. 

It is interesting to note the similarity between SOC 
models and nonequilibrium lattice automata with multi- 
ple adsorbing states j40| . These models present a critical 
phase transition separating two regimes: above the tran- 
sition there is a finite density of active sites, while below 
the transition point this density is zero and the system 
freezes in one of the many stable configurations. In the 
following, using the formalism developed for this class of 
models we will make this analogy more precise. 

III. DYNAMIC MEAN FIELD APPROXIMATION 

In the previous section we generalized the definition of 
SOC automata to the fast-driving regime, thus removing 
the assumption of timescale separation commonly em- 
ployed in simulations. This will turn out to be particu- 
larly convenient since the restored locality of the dynam- 
ical rules allows a simpler description of the models. 

The most generic description of SOC models is through 
a d-dimensional stochastic cellular automaton with N — 
L'^ sites, where L is the lattice size. Each site i on the lat- 
tice is characterized by an occupation variable ai which 
can assume q different values: for instance the possible 
energy levels Zi or the three different FF model states. 
The complete set a — {ai} of lattice variables specifies 
a configuration of the system. The dynamical evolution 
of the system is determined by the transition probability 
W{(J I a') from the configuration a' to the configuration 
a. At each time step the state of a given site depends 
only on the previous state of the site itself and the set 
of sites interacting with it. The most general transition 
probabilities in the homogeneous and symmetric case is 

N 

W{a\a')^l[w{a,\a') (4) 

i=l 

where w{ai \ a') is the one site transition probability 
that depends on driving and reaction rates. The single 
site transition probability should satisfy the normaliza- 
tion property 

Y.w{a.,\a') = l. (5) 

Because of the intrinsic nonequilibrium behavior of these 
systems, we have to consider the time-dependent prob- 
ability distribution P{a,t) to have a configuration a at 
time t. From this distribution we can compute the aver- 
age value of any function of the state A{a) 

(A(t))=^A(a)P(a,0. (6) 

a 

The time evolution of the probability distribution is gov- 
erned by the master equation (ME), which in continuous 
time reads as 
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t) = J2 W{cj I a')P{cj', t) - W{a' \ a)P{a, t). 



(7) 



The specific form of W determines the dynamics of the 
model and the steady state distribution. Typically, SOC 
systems show a stationary state in which all the single 
time averages are time independent. To this state cor- 
responds a stationary probability distribution P((t) = 
P{<J,t —>■ oo). For equilibrium systems the stationary 
distribution has the Gibbs form P{a) ~ exp(— /3i/(cr)), 
where H{a) is the Hamiltonian. For SOC systems, like 
other nonequilibrium systems, there is not such a gen- 
eral criterion, but we have to solve explicitly the ME 
in the stationary limit. In practice, this is a formidable 
task which is accomplished just in very few cases. It is 
then necessary to use approximate methods in order to 
describe the collective behavior of these systems. The 
simplest available method is the dynamic cluster varia- 
tion approach, which involves a hierarchy of evolution 
equations for the probability distribution of configura- 
tion of cluster of k sites: Pk{cri, • ■ • , Cfc). If the system 
is homogeneous, the distribution of cluster of k sites will 
be position independent. It is easy to recognize that Pi 
represents the average density of sites in a certain state, 
while Pk>i characterize the correlation properties of the 
systems. Unfortunately, the time evolution equations for 
each of the Pk depends on the higher correlation func- 
tions: the dynamical equation for the average densities 
depends on the two point correlation functions, the two 
points correlations on the three point correlations and so 
on. We have therefore an infinite chain of coupled equa- 
tions. The dynamical mean-field approximation consists 
of neglecting correlations up to a certain order. In the 
n-sites approximation cluster probabilities are decoupled 
as a product of n-sites probabilities. This approxima- 
tion has proved to be quite instructive for a qualitative 



description of the critical behavior of nonequilibrium sys- 
tems IQ. 

Before proceeding to discuss in detail the single-site 
MF approximation for stochastic SOC models, we first 
discuss the basic symmetries of these systems, which will 
play a fundamental role in formulating a common de- 
scription. We can reduce the number of states each site in 
the system can assume, noting that we can always iden- 
tify three main states: stable (ai = s), critical [ai = c) 
and active {at = a). Stable sites are those that do not 
relax (become active) if energy is added to them by exter- 
nal fields or interactions with active sites. Critical sites 
become active with addition of energy. Active sites are 
those transferring energy; they are interacting with other 
sites (usually nearest neighbors). Indeed, SOC refers al- 
ways to systems in which the only state that generates 
dynamical evolution is the active one; i.e. stable and crit- 
ical sites can change their state only because of external 
fields or by interacting with an active nearest neighbor. 
Therefore, SOC model correspond to three states CA on 
d-dimensional lattices. This description is only approx- 
imate, since a certain amount of information is lost in 
grouping together stable sites. For instance, in the BTW 
model [|l| we have several energy level which pertain to a 
stable site, but we can take this fact into account intro- 
ducing in the ME some effective parameters. The three 
states description is exact for the Manna model |^ and 
the FF model 1^. 

In the simple MF single-site approximation, we denote 
by Pa,Pc and ps the average densities of sites in the ac- 
tive, critical and stable states respectively. In the case of 
homogeneous systems, these densities can be written as: 



p„(t)=.^<5(a,-.c)P(a,t). 



(8) 



The dynamical equations for the average densities are 
obtained from the ME by using Eq. (0) 



The above equation can be simplified by using the nor- 
malization condition for the transition probabilities 

J2Y[w{a',\a) = l, (10) 

{a'} t 

Y^Siaj - K)l[w{a, I a') = w{aj = k \ o'). (11) 

{a} 

Equation (^ can be further simplified when the interac- 
tions are only among a finite set of sites. In this case, 
with o' = {<y'iT<y'ij^e} ^6 denote the site i and the set of 



sites that can interact with it — usually a finite num- 
ber of sites or more commonly just the nearest-neighbors 
(n.n.). By restricting the sum and dropping the site in- 
dex, because of the homogeneity, we finally get 

{a'} 

It is worth to remark that in the above expression the 
set a' = {cr-, cr-_|_e} refers to the generic set of interact- 
ing sites which depends upon the particular dynamical 
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rules and lattice geometry. In presence of a non-local in- 
teraction, the set a' can correspond to the entire system. 
This presents a very difficult problem that can be treated 
introducing a suitable regularization. 

In general, we have therefore that the evolution equa- 
tions of the average densities are still coupled to the prob- 
ability distribution of configurations of a set of interact- 
ing sites. In order to have a set of closed equation for the 
densities we truncate the evolution equations by using 
a dynamical MF single-site approximation. In this MF 
scheme we approximate the probability of each configura- 
tion a as the product measure of single-site probabilities 

i i 

thus neglecting all correlations in -P(cr). Introducing this 
approximation in Eq. (^, we obtain the MF reaction rate 
equations which depends just on the single-site densities 
and can be symbolically written as 

d 

-^Pn ^ F,,{pa,Pc,Ps) K=a,C,S, (14) 

where depends upon driving fields and interactions 
parameters through the transition rates w. In addition, 
because the densities must preserve normalization, two 
of the above equations supplemented with the condition 
Pa + Pc + Ps = Ij are enough to describe completely the 
system. 

In practice the form of the rate equations depends upon 
the specific model. Nevertheless, we can write the gen- 
eral structure of the equations describing SOC models by 
simple considerations. In general, Ff^ can be expanded as 
a series of the average densities: 

n n,£ 

where the constant term is set to zero in order to get 
a stationary state. The first order terms are the tran- 
sition rates generated by the external driving fields or 
by spontaneous transitions. The second and higher or- 
der terms characterize transitions due to the interaction 
between different sites. In SOC models, only the active 
state generates a non trivial dynamical evolution, while 
stable or critical sites can change their state only be- 
cause of the external field or the presence of an active 
n.n. site. Since the critical point is identified by Pa = 0, 
in correspondence with a vanishing external field, we can 
neglect second order terms in the density of active sites. 
The solutions of the stationary equations (^/Ok =0) are 
function of the effective parameters /",/"'^, which de- 
pend on the details of the model. It is expected that the 
critical behavior is not affected by the specific values of 
the parameters, while universality classes will depend on 
constraints imposed on the equations because of symme- 
tries and conservation laws. 



IV. MEAN-FIELD ANALYSIS OF SANDPILE 
MODELS 

We consider here the explicit application of the single 
site MF approximation to the class of sandpile models. 
A simpler derivation based on symmetry considerations 
can be found in Ref. jl^ . In the previous section we have 
shown that the MF dynamical equations reduce in this 
approximations to the the following expression 

§iPnit) = E "'(^ I n p< (^) - '^-(^)- (16) 

{a'} 

where a' denotes the set formed by a single site and its 
set of interacting sites as specified by the dynamical rules. 
All the dynamical information of the system is contained 
in the transition rates w{k | ct'). Unfortunately, the sand- 
pile model is inherently non-local because of the implicit 
timescale separation. A site can receive energy only if the 
system is quescient. This implies that transition rates 
depend upon the whole set of lattice variables present in 
the system, giving rise to a strongly non-local dynamical 
rule. In order to treat the model analytically we have to 
regularize this interaction by a suitable parametrization 
which allows to recover the non-locality in some particu- 
lar limit. The simplest regularization has been discussed 
in Sec. II, introducing the external flow of energy added 
to the system. We describe this external flow by the 
probability per unit time h for a site to receive a grain 
of energy. The transition rates are now local, depending 
only on the field h and the state of the nearest neighbor 
sites ((T„.„.) that determine the toppling dynamics. The 
total amount of energy added to the system at each time 
step will be J„ = hL'^ ||. The non-locality of the dy- 
namical rules is recovered in the limit h (see Sec. 
VI), that corresponds to an infinite timescale separation. 
The external field h has been historically introduced in 
Ref. . Unfortunately, these early papers failed to ad- 
dress consistently the role played by driving and conser- 
vation, and led to several inconsistencies (see Sec. IV B). 
Other regularization can be introduced as well in the ME 
treatment. We limit our discussion to the present one for 
reasons of simplicity. Nevertheless a more accurate char- 
acterization of the degree of non-locality actually present 
in the infinitely slowly driven sandpiles can be obtained 
via more refined regularization schemes |^6|] . 

Since locality is restored {a' = {^i, o-^.„.}), we can de- 
rive the MF equations for the density of active sites by 
considering the leading order in h and pa in Eq. (|l^). 
The transition rates obey u;(a | a,an.n) — 0, because 
an active site always transfers its energy, thus becoming 
stable at the next time step independently from its n.n. 
sites. In this way, we are neglecting higher contribution 
due to the presence of multiple active n.n sites, which can 
transfer energy to the active site sustaining its activity. 
The only allowed transition to the active state are due 
to critical sites which receive energy from the external 
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driving or from active n.n. sites. In the absence of active 
n.n sites, we have w{a \ c, an.n. ^ a) = h. We can then 
obtain the contribution to the dynamical MF equation 

X! I C, tr^.n. 7^ W Pa'^{t) ^ hpc{l - Pa)^ , 

Wn n } i£n.n. 

(17) 

where Z represent the lattice coordination number; i.e. 
the number of n.n. sites. If the sites does not receive 
energy from outside we have to consider the possibilities 
that one of the n.n. sites is active and transfers energy 
to it. This process corresponds to 

w{a \ c,(Ti= a, aj^i ^ a) = [l - h)^{l - p), (18) 

where i,j(z n.n.. The right term represent the probabil- 
ity that a critical site receive an energy grain only from an 
active n.n.. This is equal to the ratio between the num- 
ber of sites g involved in the dynamical relaxation process 
and the total number of n.n.. For instance, g — 2d for the 
BTW model ^ or g = 2 for the Manna model [|9|. In 
addition, we have to consider the probability p that the 
active site does not transfer its energy because of intrinsic 
dissipation or because it is a boundary site. The above 
transition rate is valid only for homogeneous processes 
and therefore excludes directed models. The total con- 
tribution due to this process, considering the multiplicity 
of active n.n. sites, is given by 

^ w{a I c, cr,' = a, cr^-_^, ^ a)pcPa J| Pa'^ {t) 

^{g^e)p,pa{l^h){l-paf'\ (19) 

where the parameter e conveniently identifies the aver- 
age energy dissipated gp in each elementary process. It 
is worth to remark that e is present also for fully conserva- 
tive systems, being an effective term due to the boundary 
dissipation: it acts as an external tunable parameter in 
the case of bulk dissipation and accounts for size effects 
in finite systems. 

Neglecting higher orders in h and pa from the Eqs. ( |l7| ) 
and (|l9|) , we can finally write the MF dynamical equation 
for the densities of active sites 

d 

J^Pa{t) = -Pa{t) + hprXt) + 

{g^e)p,{t)pa{t) + 0{hpa,pl). (20) 

Next, we derive the dynamical MF equation for the 
density of stable sites, following the same strategy used 
above. Since, at lowest order, active sites become sta- 
ble with unitary rate, we have that w{s \ a,a'^ = 
1 -|- 0{hpa, Pa), yielding a contribution to the MF equa- 
tion which is 

^(■S I ^'^n. nJPa (*) = Pa + C(/ipa , Pa ) • 

{^'n n} i€n.n. 

(21) 



Since critical sites never become stable, we have also that 
w{s I c,<„ ) = 0. 

Energy conservation imposes that energy is stored in 
stable sites until they become critical. This implies a 
non unitary rate w{s \ s,a'^^ ). The simplest way to 
derive this term make use of the normalization condi- 
tion that yields w{s \ s,cr^ „ ) = 1 — w{c \ s,a'„ „ ). 
In fact, these transition rates are nearly equivalent to 
those from critical to active sites. The only difference 
is that only a fraction u of stable sites receiving an en- 
ergy quantum will contribute to the s — s- c process, i.e. 
the fraction of stable sites which are sub-critical. There- 
fore the reaction rates are related by the factor u as 
w{c I s, a'„ „ ) = uw^a \ c, cr^ „ ). Recalling the deriva- 
tion of Eq. (E3) , it is straightforward to obtain 

{l-w{c\s,a'„^J)ps Yl P<yfi) = 
Ps - uhps - u{g - e)psPa + 0{hpa,pl). (22) 

Adding all these contributions, wc finally obtain the dy- 
namical MF equation 

d 

— Ps(t) = pa{t) - Uhps{t) 
Ot 

+u{g - e)ps{t)pa{t) + 0{hpa,pl), (23) 

that together with Eq.(|20|) and supplemented with the 
normalization condition fully describes the MF evolution 
of sandpile automata. 

In deriving the MF equations we have made an ap- 
proximation, introducing the parameter u to take into 
account the presence of several energy levels instead of 
a single stable level. For three levels models u = 1 and 
this description is exact, while for multilevel models like 
the BTW ^ the parameter u can be determined self- 
consistentlyin the stationary state using the energy con- 
servation |^3[. Here we show that we can also obtain u 
by the full description of the MF equations. We consider 
a generic sandpile model in which the energy threshold 
is Zc and, after a relaxation event, g energy grains are 
transfered to randomly chosen neighbors. For instance, 
the d-dimensional BTW has Zc — 2d and g = 2d, but 
we can think to arbitrary values for these parameters. 
We can describe in more details these systems by intro- 
ducing the densities p„, describing the probability that 
a site is in the level n. We then have that pc = Pz^-i 
and ps = X]ra=o P"^' dynamical evolution can be 

simplified noting that an active site with energy Zc be- 
comes stable and its energy becomes n = Zc — g. One 
can show that stable sites with energy levels lower than 
n, have a zero stationary density. Without loss of gener- 
ality, we can therefore assume that the zero energy level 
is n = Zc — g . By rescaling the energy levels in this way. 
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we get pc = Pg-i and Ps = J2n=oP^- ^^'^ intermediate 
level are described by the following set of MF equation 



hL" 



Jout — ^PaL'^- 



(29) 



d_ 
di 



Pn{t) = -hpn(t) - {g - e)pn{t)pa{t) + hpn-l{t) 

+ {9 - e)pn-lit)pait) + Oihpa,pl), (24) 



where 1 < n < g — 2. In the stationary state we 
obtain pn = Pn-i ■ ■ ■ = Po and, noting that u ~ 
Ps-2/(Z]n=CL£n)' recover the result u = 1/(5 - 1) 
obtained in |13|. This result expresses the energy conser- 
vation and fixes the stationary solution consistently with 
the energy balance. Noticeably, in the Manna model ||49| ] 
for which 5 = 2, we obtain u = 1 as it must be for a 
three states model. As a last remark we point out that 
summing up the above set of equations with the one for 



d 

OlPoit) = -hpo{t) + 

{9 - '^)Po{t)Pa{t) + Pa{t), 



(25) 



we obtain Eq. ( pS]) as a function of the parameter u. 

When the system is far from the stationary state, the 
parameter u will be in general time dependent. In the 
following we will always consider stationary properties 
or homogeneous perturbations which leave u unchanged, 
but we could think of situations in which u has not its 
stationary value. This corresponds to a systems kept far 
from its "natural" configuration. This can have strong 
influence even on the critical properties of the system as 
in the case of CA with absorbing states. These features 
will be discussed elsewhere p6| . 

To study the stationary MF solutions, we consider 
the simple three level case and we determine u self- 
consistently as in Ref . . Combining the normalization 
equation with the stationarity limit of Eqs. ( pO| ) and |2^, 
we obtain the set of equations 



Pa = hpc + [g - e)pcPa: 
Pa = uhps + u{g - e)psPa, 

Pa = 1 - Ps - Pc- 



(26) 



After some algebra from Eqs. (|2^) we obtain a closed 
equation for pa 

u{g - e)pl + {l + u{l + h~g + e))pa -~uh = 0. (27) 

We can expand Pa{h) for small values of the field h. The 
zero order term in the expansion vanishes and we obtain 
a leading linear term: 



Paih) 



uh 



l + u{l-g + e)' 



(28) 



This result has to be consistent with the global conser- 
vation law, which states that the average input energy 
flux Jin must balance the dissipated flux Jout- In the 
stationary state the conservation law can be written as 



By comparing Eq. (g8|) with Eq. (^9|) we obtain that 
u — 1/(5 — 1), which is the result we have previously 
obtained from the complete analysis. In the limit /i — > 
the densities are therefore given by 



h 



Pa 



0{h), 



1 
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0{h). (30) 



The numerical values for the density of critical and stable 
sites are non-universal quantities and depend on the lat- 
tice geometry and dynamical rules of each speciflc model 
via the parameter g. The result for pc can be directly 
compared with the estimates from numerical simulations 
of several models by substituting the correct value of g. 
For the original BTW model {g — 2d), extensive numer- 
ical simulations on the density of energy levels can be 
found in Ref |3^. As expected the agreement with the 
MF result increases for high dimensional systems and we 
recover the exact result in the limit d — > 00. 

We next discuss the critical behavior of these systems. 
The balance between conservation laws and the dissipa- 
tion are essential for the critical behavior of the model, 
as also pointed out in The model is critical just 

in the double limit h,e 0, h/e — > 0, similarly to the 
forest-fire model We are going to see that in this 

limit the zero field susceptibility of the system is singu- 
lar, signaling a long-ranged (critical) response function. 
The onset of the critical behavior is then recovered in 
the limit of vanishing driving field corresponding to the 
locality-breaking in the dynamical evolution. In analogy 
with non equilibrium phenomena p0| , p2| , the one par- 
ticle density of active sites is the order parameter and 
goes to zero at the critical point. The driving and dis- 
sipation rates identifies the two control parameters, i.e., 
the relevant scaling fields. We can then distinguish dif- 
ferent regimes as a function of the control parameters. 
The model is supercritical for h > and e > h, while 
for /i — > and e > it is subcritical and the dynamics 
displays avalanches. The phase diagram is somehow sim- 
ilar to that of usual continuous phase transitions, if we 
replace h by the magnetic field and e by the reduced tem- 
perature (see Fig. We can fully exploit this analogy 
by allowing the parameter e to assume also negative val- 
ues. This correspond to a sandpile in which a positive net 
amount of energy enters the system during the avalanche 
activity (negative dissipation). The resulting supercriti- 
cal regime is analogous to many nonequilibrium systems 
with negative reduced control parameter. However, for 
h > e the system has only trivial stationary state, since 
Pa would have to be greater than one to satisfy Eq. (|2S|). 
Thus, in the non-trivial supercritical region h and e can 
not be varied independently because the global conserva- 
tion imposes that h < e. This restricts the scaling behav- 
ior to particular limit values of the control parameters. 
In the following we individuate the regimes correspond- 
ing to the standard sandpile numerical simulations and a 
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completely new scaling regime in the supercritical region 
of the phase space. 



A. The subcritical regime 



Hence, at large r we have the solution y(r) ^ r^""*. A 
similar result has been obtained in Ref . |20[| . In the pres- 
ence of boundary or intrinsic dissipation, the system ac- 
quires a finite correlation length and we can establish the 
general scaling form 



The standard numerical simulations of sandpile mod- 
els are carried out in the presence of an infinite timescale 
separation. As already discussed in the previous sections, 
this implies an infinite slow driving of the systems, i.e. 
/i — > 0. In this limit there is a single control parameter, 
the intrinsic or boundary dissipation e, and the order pa- 
rameter Pa is identically zero in the steady state. To 
describe quantitatively the critical behavior, we study 
the effect of a small perturbation Ah on the steady state 
density 

Apa{x,t) ^ Xh.Ax - x';t -t')Ah{x',t')dx'dt' (31) 



where Xh,e{x — x';t — t') is the response function of the 
system. We define the total susceptibility Xh.n of the sys- 
tem as the integral over space and time of the response 
function, and it is shown in appendix A that for the sta- 
tionary state 



lim Xh,. 



dpa{h) 



dh 



h=0 



This immediately gives the zero field susceptibility 



(32) 



(33) 



which diverges as e ^ 0. The system is in a subcritical 
state for any value of e different from zero. For e = 
the system reaches a critical point in which the response 
function becomes long-ranged and the susceptibility di- 
verges. Close to this critical point the scaling behavior is 
characterized by the scaling laws Xe ~ ^"''^i with 7 = 1, 
and by the divergence of the correlation length ~ e""^. 

An important result can be derived from the response 
function by defining 



Xe{r) = / Xe{r,t)dt 



(34) 



as the average total response received at position r, 
when Ah{x',t') = S{x') Since energy is transfered lo- 
cally and isotropically, the net energy current is given by 
j ~ dx{r)/dr. For locally conservative models the en- 
ergy current j must satisfy in average the conservation 
law 



jda = cost, 



(35) 



where da is the d — 1-dimensional surface element. This 
ensures that the energy flowing into the system is bal- 
anced by the dissipated energy in the stationary state. 



Xeir) = r(r/C), 



(36) 



where T(r/^) is a cut-off function for r >> ^. This im- 
mediately gives the following relation between ^ and the 
zero field susceptibility 



Xe 



= / xe{ry-'dr ~ e 



(37) 

= 1/2 



We find the MF value of the correlation exponent 1/ 
by substituting ^ ~ and comparing with eq.(|33|). 

We can use these exponents to characterize the finite 
size scaling of the conservative sandpile model, since our 
MF analysis treats both boundary and bulk dissipation 
in the same way. In conservative systems, when the size is 
increased the effective dissipation depends on the system 
size and we assume that e ~ L~'^. In fact, the dissi- 
pation rate is given by the probability to find a border 
site instead of a bulk site during an avalanche. Thus, the 
exponent links the dissipation rate with the system 
finite size, providing a unified view of locally dissipative 
and open-boundary models. In the conservative case, the 
characteristic length of avalanches should go like £, ^ L 
to ensure dissipation of energy outside the boundaries. 
This implies the scaling relation = 1, which imme- 
diately gives (1 = 2. We show in appendix A that the 
susceptibility scales as the average avalanche size 



which implies x 
scaling law 



Xe s >, (38) 
2^M7 From this result we obtain the 



(s) ^ for L ^ 00, 



(39) 



which has been found numerically in various dimensions 
p6[ and and rigorously proven in d = 2 ||T^. Our 
explanation implies that the diffusive behavior of the 
avalanches is due to the global conservation law. 

Summarizing all these results we obtain a first set of 
MF exponents 



7 = 1, (1 = 2, iy = l/2. 



(40) 



In deriving these exponents we made only use of conser- 
vation laws, therefore we expect that these results hold in 
all dimensions. In the next section we will confirm these 
results by numerical simulations performed in the fast 
driving regime and we will discuss some results already 
published in the literature. 

In the subcritical regime the dynamics takes place in 
the form of avalanches, but if ft. = the system rapidly 
decays in one of the adsorbing configurations; the ones 
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with no active sites. All of them are stable in absence of 
the driving field. It is useful to characterize the prolifer- 
ation of active sites starting from a seed initial condition. 
In close analogy with CA with adsorbing states, we study 
the spreading of active sites after a small perturbation. 
We prepare the system in an initial state consisting of 
a single active site, i.e. an infinitesimal perturbation in 
the driving field Ah{x,t) = S{t)S{x). Since h{t > 0) = 0, 
active sites caimot be produced spontaneously from crit- 
ical sites and can only appear due to the spreading of the 
initial perturbation. The properties of this process close 
to the critical point characterize the avalanche behavior 
typical of SOC phenomena. 

Following Grassberger and de la Torre we consider 
the probability that a small perturbation activates s sites 
(an avalanche in the SOC terminology) 

P{s,e)^s-^g{s/sc{e)), (41) 

where Sc ~ e"^/'^ is the cutoff in the avalanche size. The 
perturbation decays in the stationary subcritical state as 

Pa{t)^t'^T{t/t,{e)). (42) 

Here tc denotes the characteristic time which scales as 
tc - e~^. We can also introduce the scaling exponents 



which relates cutoff lengths to the characteristic size, 
Sc ~ and tc ~ ■ These exponents are related by 
the following scaling laws 

D = —, zi' = A. (43) 

1^(7 

Another scaling relation between critical exponents can 
be obtained from Eq. ( ^ ) (see Appendix A) 

Xe s J s-^+^g{slsc{e))ds ~ e(--2)/-, (44) 

which implies 

7 ^ (45) 

cr 

To obtain the MF values of the avalanche exponents, 
we solve the evolution equation for a small perturba- 
tion close to the stationary state. We consider pK.[t) — 
Pk + Sp^it), where Spi^{t) is the deviation of the densities 
from their stationary value. By considering small pertur- 
bation around the stationary state, keeping only linear 
term in 5pi^{t)^ and using the normalization condition we 
obtain 



— (5/9a(i) = -Spa{t) + hSpc{t) + {g ~ e)pc5pa{t) + {g ~ e)paSpc{t) 
ot 

d h € €. 

— Sps(t) = +Spa{t) -5ps{t) H -psSpa{t) H -paSps{t) 

Ot g-l g-l .g-1 

5pa{t) = ~5pc{t)-5ps{t). (46) 



In the subcritical regimes {h 0) we only keep in 
these equations the leading terms in e. Substituting in 
Eqs. (|4^) the densities given by the solution of the sta- 
tionary equation for /i ^ (i.e. Pa = and pc — I/3), 
we finally obtain the evolution equation in diagonal form 



d e 

— (5pa(i) = 5pa{t), 

dt g 



The solution of Eq. (47) is given by 



6pa{t) ^ cxp{~et/g) 



(47) 



(48) 



which implies rj — 0. The last equation also defines the 
characteristic relaxation time for an infinitesimal pertur- 
bation to be tc = g/e, yielding A = 1. We compute 
the remaining exponents using a further scaling relation, 
which we derive in appendix B, 



(r-1) 



(49) 



It is worth to remark that Eq. (^^ is valid only in 
MF theory. By combining these relations with those of 
Eq. (^ , we get the second set of MF critical exponents: 



D = 4, T = 3/2, cr = l/2. 



(50) 



It is worthwhile to remark that the numerical value of 
these exponents is the same as in other MF approaches, 
but their significance is completely different, being de- 
fined with respect to a different scaling field. All sandpile 
models with the same dynamical MF equations share the 
same critical exponents and belong to the same univer- 
sality class. However, the degree of universality is highly 
overstated, as usually happens in MF approaches. In 
particular, the exponents do not depend on the dimen- 
sionality d. The exponents describing avalanche distri- 
butions in low dimensional systems, in general, will not 
agree with the results of MF theory. For instance, it is 
still controversial if the BTW and the two-state model 
are in the same universality class. While it was believed 
for some time that this was indeed the case, recently 
large scale numerical simulations questioned this state- 
ment iTil. 
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To compute the value of critical exponents below the 
upper critical dimension, we have to use renormalization 
group techniques, which allow the correct treatment of 
the scale- free fluctuations present at the critical point. 
The renormalization group approach to non-equilibrium 
system presents several difficulties which can be overcame 
by suitable approximations ]l^ , p^ . 

B. The supercritical regime 

The supercritical region is characterized by a finite 
density of active sites; i.e. a non-zero order parameter. 
Close to the critical point, the supercritical region corre- 
spond to the parameter range h << 1, e << 1 and h ^ e. 
In this regime the order parameter is linear in h 

Pa - h^": 5 = 1, (51) 

as we obtain from Eq. ( |3C|). The same result has been 
also conjectured in Rcf. |}46{. This is analogous to the 
MF results obtained for contact processes and other non 
equilibrium CA [^o|-^, but it is in contrast with pre- 
vious MF approaches for sandpile models |^,^, which 
yielded 6 = 2. This latter incorrect result is due to an 
inconsistency present in those studies. The scaling is 
expressed in terms of the average energy (z) = J2i Pi^i 
which is treated as an independent control parameter. 
As we have just shown, (z) and h are not independent 
in the stationary state. The stationary probability dis- 
tribution of heights pi is indeed a function of the driving 
rate. Moreover, (z) can not be considered as the con- 
trol parameter even for h = 0, since it does not deter- 
mine completely the state of the system: the same value 
of (z) describes several states corresponding to different 
values of densities pi. This is a typical property of CA 
with multiple absorbing states |12[ . In preparing an ini- 
tial condition consisting of a localized active region, one 
has considerable freedom to choose the initial state (the 
adsorbing configuration). In order to observe the criti- 
cal properties of the system, we should chose one of the 
"natural" initial conditions, which can be obtained by 
the dynamical evolution for infinitesimal driving in the 
long time limit. In numerical simulations this is equiv- 
alent to first prepare the system in the stationary state 
in presence of the timescale separation, and then average 
over the many different realization of the avalanches. 

The exponent S has been dcffncd in previous works in 
analogy with usual continuous phase transitions, where it 
characterizes the scaling of the order parameter in pres- 
ence of an external field when the other critical param- 
eters are set to zero. In SOC, however, h and e are not 
fully independent because the of the e > h condition. 
Thus in the supercritical regime, the scaling of physical 
quantities is a general homogeneous function of the fol- 
lowing kind 

f{h,e) = b^f{by-h,by'e), (52) 



where 6 is a scaling factor and e can not be set to zero 
for finite h. The scaling behavior of the system is there- 
fore particularly complex, as in the case of conventional 
phase transitions in the region where two control parame- 
ters are different from their critical values. We can define 
a consistent scaling regime with respect to the reduced 
variable = h/e in the double limit h and e — > 
with the supplementary conditions that e <<((><< I. 
These limits define a parameter region which is identified 
roughly as e << 1 and « h « e. In this region, 
the MF approximations shows that the order parameter 
is positive and scales as 

Pa = ^^ (53) 

with (3=1. The exponent (3 characterizes the scaling 
behavior of the order parameter with respect to the re- 
duced parameter </> and should not be confused with the 
exponent S. In order to uncover the scaling of the char- 
acteristic lengths of the system with respect to the pa- 
rameter 0, we study the evolution of small perturbations 
around the stationary state. As in the previous section, 
we denote by Sp^it) the deviations from the stationary 
state and write the dynamical equation keeping only the 
leading terms. We thus neglect in Eqs. (|6|) the terms in 
h and e, and keep the terms in (j). For this we compute 
the first order correction in (j) to the values of the station- 
ary densities. These results to be pa = 4>, Pc = — 4')/ 9 
and Ps = {g — 1)(1 — 4')/9i ^^'^ by substituting in the 
dynamical MF equations we get 

d 

— (5pQ(t) = -{g - l)(j)Spa{t) - g(j)Sps{t) 

6psit) = +4Spait) - ^cj,6ps{t) (54) 
dt g-1 

By diagonalizing Eqs. (p^, we find the eigenvalues 

A± = -0/±(.g) (55) 

where f±{g) has always positive real part. Both eigen- 
values are thus negative and linear in (j) and represent 
the inverse of the relaxation timescale of a perturba- 
tions around the stationary state. We have therefore that 
tc ^ 4>~^ , with A' = 1, implying that the characteristic 
time scales as in the subcritical regime. The solution for 
the spreading of the density perturbation has the form 
Spa(t) ~ !F{t/tc), yielding as in the subcritical regime 
r]' = 0. 

Eqs. ( [54| ) describe how a locahzed perturbation decays 
in the stationary state. As in the previous section, this 
decay can be related to the susceptibility 

X^- j 5pa{t)dt <i>-^' (56) 

with 7' = 1. A characteristic length ^ is associated 
with the characteristic time of ffuctuations. Since en- 
ergy is transferred homogeneously and isotropically, we 
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have that ~ C^, as in the subcritical regime. By com- 
paring this relation with Eq. (|5^), we obtain ~ 0^^, 
or ^ (f)^'^ with v' = 1/2. We can obtain a clearer 
picture of this behavior using the avalanche representa- 
tion. The condition that the time between two energy 
addiction (the driving timescale) is much longer than the 
fluctuation timescale can be written as 

h « - 4>. (57) 

This condition is implicitly verified in the limit we are 
considering [h << e << cj)). Under this assumption it 
is very unlikely for fluctuations to overlap. Thus, on av- 
erage, each event is separated from the others and can 
be defined as in the subcritical case as an avalanche. In 
the stationary state the average non-zero order param- 
eter is produced by the random appearance of a finite 
number of non-overlapping avalanches. Furthermore, we 
can identify the response function with the time evolu- 
tion of an avalanche, and we recover in the limit /i — > 
that x<t> ~< s >■ On its turn the latter implies that in 
this regime the average sizes of avalanches diverges as 
< s >~ as — > 0. 

These results are particularly interesting because it is 
the first time that the supercritical scaling behavior is 
characterized by the parameter (j) = h/ e. The MF predic- 
tions, also if not exact from a quantitative point of view, 
should refiect the right symmetry of the problem. Thus 
we expect that a supercritical regime should be found 
also from numerical simulations, but with different val- 
ues of the critical exponents. Although the analytical 
formulation of the supercritical sandpile is quite simple, 
the actual implementation of the numerical simulations 
is not straightforward since the scaling regime is obtained 
only in the double limit /i — > and e ^ 0. It is not possi- 
ble to use anymore the timescale separation, because we 
want investigate the behavior for finite (p, when a single 
timescale represents the system. Work is in progress to 
obtain a numerical confirmation of the MF predictions 
for the supercritical regime. 

Usually the equilibrium version of the fluctuation- 
dissipation theorem allows the calculation of the upper 
critical dimension of models by requiring the consistency 
of the MF theory. The theorem links the susceptibilities 
with the equilibrium fluctuations from which Ginzburg- 
like self-concistency criteria are derived. The existence of 
similar theorems for nonequilibrium steady-states would 
be clearly useful. From general considerations, it is pos- 
sible to show that fluctuation-dissipation-like theorems 
do exist for nonequilibrium models. However, the phys- 
ical quantities that satisfy these theorems are generally 
unknown, and their construction requires the complete 
solution of the system's dynamics. This fact is very 
unsatisfactory from a practical point of view and has 
given rise the terminology of "violation of the fluctuation- 
dissipation theorem" . For stochastic SOC we are there- 
fore unable to use the susceptibility in order to estimate 
the stationary fluctuations. This does not allow, at the 



present level of approximation, to find the upper critical 
dimension dc of this class of models. From the theoretical 
point of view the situation is still very controversial. Sev- 
eral theoretical estimates give dc = 4 , that has been 
obtained also from recent numerical simulations p8[ . In 
contrast, other numerical studies and the similarity 
to percolation led several authors to conjecture — Q. 
Despite the fact that exponents have the same numerical 
values, our approach points out that the MF theory for 
SOC models is very different from that of percolation. 
For instance, we have shown that the exponent 5 is not 
well defined in our case. We can understand therefore 
that some scaling relations such as 5 = (r — 2)~^ which 
are valid for percolation do not apply here. 



C. Numerical simulations 

In this section we compare the results of numerical sim- 
ulations with the prediction of MF theory, and in partic- 
ular we expect that the set of MF exponents related to 
the global conservative nature of the model are valid also 
in the low dimensional cases. Sandpile models have been 
extensively studied only in the subcritical state [|5|-Q. 
Most of the numerical results refer to avalanche distribu- 
tion in the conservative limit, with open boundary condi- 
tions. In these conditions the finite size scaling has been 
found to be problematic and despite the use of very large 
scale simulations there is not a complete agreement on 
the values of the exponents |^^. The reason for this is 
probably that open boundary conditions impose a value 
for the effective dissipation which depends on the lattice 
size and does not act homogeneously through the system. 

We simulate numerically the BTW model with finite 
driving rate h and boundary dissipation in rf = 2. In 
this case the dissipation is implicitly considered through 
the open boundary conditions. When a boundary site 
topples, it dissipate part of the energy outside, without 
transferring it to the neighbors. The driving rate h is 
introduced as the probability for unit time that a site re- 
ceives an energy grain. Apart from the driving the simu- 
lation proceeds as in Ref. Q. We see in Fig. || that the 
density of critical sites goes to zero linearly with h{5 — 1) 
with a slope that increases with the system size as . 
This is in agreement with the MF theory which predicts 
that the susceptibility scales as L'^''', with = 2. 

To observe more clearly the scaling with dissipation of 
the sandpile model we study the BTW model with pe- 
riodic boundary conditions and fixed dissipation e. We 
model the dissipation introducing a probability p = e/g, 
for which the energy in a relaxation event is lost, instead 
of being transfered. In Fig. || we plot the susceptibil- 
ity Xe = dp/dh as a function of e. We observe the 1/e 
behavior (7 = 1) predicted by mean-field theory. We 
should add to this the value of = 0.5 that was obtained 
studying a dissipative sandpile model in d = 2 [5C| ] . 

In summary, we have shown that some MF features 
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are present in the d = 2 sandpile model. The global 
conservation law imposes the exponents /i, v and 7 to as- 
sume their MF values. This strongly support the MF 
picture provided here. The other critical exponents do 
not necessarily assume mean-field values in d — 2 and 
we do not analyze them here. Extensive measures of 
these exponents can be found in the literature ps]-^. 
We are currently performing simulations of the model in 
the supercritical regime in d > 2 and the results will be 
published elsewhere Q. 

We expect the critical properties, such as exponents 
and scaling, to be universal for the class of "natural" 
configurations. It is interesting to note that spreading 
exponents could be sensitive to initial configurations dif- 
ferent from the "natural" ones. This has been extensively 
studied in noncquilibrium phase transitions, where the 
critical behavior is changing with respect to the adsorb- 
ing configurations used as initial conditions. In equilib- 
rium critical phenomena, the appearance of continuously 
variable exponents is usually associated with a marginal 
parameter, that can be the case for critical phenomena at 
adsorbing states. Unfortunately, no numerical data are 
available on sandpile models to place these speculations 
on a firmer basis. 



D. On the role of conservation 



in d = 2 for a > Qfc the system is critical. There is still 
not an agreement on the precise value of ac, which was 
first estimated as ac — 0.05 and later found to be 
higher (etc — 0.18) while it was claimed in Ref. |6|] 
that ac — 0. 

The mean-field analysis of this model is not easy be- 
cause of the continuous number of levels a site can as- 
sume. Using an approximate analysis of the random- 
neighbor model it was claimed in that ac — 0.22, 
a value that was found in agreement with simulations. 
A complete analysis of the master equation later re- 
vealed that ac — 0.25 (the conservative case), showing 
also the presence of very strong finite size correction [ p4| . 
From this analysis, it appears that the random neighbor 
model behaves like the BTW model. Criticality is only 
reached in the conservative case, in the limit of zero driv- 
ing rate. The situation in two dimensions is still contro- 
versial and it is believed that the inhomogeneity created 
by the open boundary conditions is responsible for the 
observed power law distributions. 

The role of conservation for criticality remains still 
open in this models, while it is now agreed that in MF 
theory conservation is a necessary condition to achieve 
criticality in sandpile models. In the next section we will 
discuss a model that shows criticality without con- 
servation even in MF theory. The price to achieve this 
result will be an additional driving rate. 



In the previous discussion, we emphasize the role of 
conservation in the dynamics of sandpile models. Using 
conservation laws, we found a subset of critical exponents 
that retain their mean-field values even in low dimen- 
sions. Conservation has an important effect on critical- 
ity as well, since the amount of dissipation plays the role 
of control parameter. Only by tuning this parameter to 
zero the system is critical. 

The role of conservation in SOC models has been the 
object of a long controversy. It was first claimed that con- 
servation was a necessary conditions for criticality in this 
class of models | |2l] , p^ . The result on dissipative sand- 
pile models seemed to confirm this conclusion . Later 
on, simulation of an earthquake model contradicted this 
results |5^. The model studied in ||5^ is a continuous 
height sandpile model in which energy in dropped uni- 
formly over all the lattice at infinitesimal rate. This in 
practice corresponds to rise the heights of all the sites by 
the quantity needed for the higher site to become unsta- 
ble. When a site i is unstable {zi > Zc) the relaxation 
rules are 

^ (58) 

Zj Zj + azi (59) 

where j are the nearest neighbors of i. The dissipation 
parameter a Q can be tuned: for a — l/(2d) the system 
is conservative. It has been observed in simulations that 



V. NON CONSERVATION AND CRITICALITY: 
THE FOREST-FIRE MODEL 

We have discussed that conservation in sandpilcs is 
crucial to achieve criticality. The controversial issue of 
continuously driven models raises the question of the pos- 
sibility that timescale separation alone can produce scale 
invariance in systems without conservation laws. In this 
context the forest-fire model acquires a very important 
role in that it is a nonconservative automaton displaying 
criticality. 

As outlined in Sec. Ill we can describe the model in 
the same language used for sandpile models We iden- 
tify burning sites with the active sites since they inter- 
act with other sites independently of the driving fields. 
Furthermore, their density vanishes in the limit of small 
driving field /. In the same way, trees correspond with 
critical sites and empty sites with stable sites. In this 
case the general three state description is exact. Using 
this language, we can emphasize differences and analo- 
gies between FF and sandpile models. While the main 
dynamical transitions are very similar, we can immedi- 
ately recognize the effect of non conservation. In the 
FF model energy in not stored and new critical sites are 
created by a second independent field, the tree growth 
probability p. Thus in the FF model we replace h ^ f 
and uh — > p. The last substitution introduces a new in- 
dependent field (or a new timescale p~^) related to the 
injection of critical sites in the system. Since energy is 
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not accumulated, there is no need of an additional dis- 
sipation, so that in the FF model there is no parameter 
playing the role of e. 

The FF model was originally introduced in the limit 
/ and p ^ by Bak, Chen and Tang in Ref. 
The model was claimed to show SOC, but later Grass- 
berger and Kantz 1 53 1 showed that in d = 2 the model was 
critical in a trivial sense. The system shows a diverging 
characteristic length that is essentially the distance be- 
tween straight fire fronts. This implies that the dynamics 
is governed by the average tree density over larger and 
larger regions. In higher dimensions the possibility of 
a non trivial behavior has not been ruled out as recent 
work seems to suggest Q. Drossel and Schwabl ||5^] 
introduced the ignition or lightning probability /. This 
field sustain fires and the system flows in a stationary 
state which shows critical properties in the double limit 
/ << p << 1. This version of the model has been the 
subject of several studies both analytical |]l2| , ^ , |65| -|6^] 
and numerical p8[-f70t. 

Despite the various efforts, the two versions of the 
model were always studied as very different cases, al- 
most two different models. For this reason, it is difhcult 
to find in the literature a precise connection among the 
two different regimes. In this section we recover within 
our framework many results already present in the liter- 
ature. By recasting these results in the language devel- 
oped for sandpile automata, we provide a unified picture 
of both models. We discuss the FF model in terms of 
the response function singularities and we show that the 
SOC-FF and the deterministic FF correspond to the su- 
percritical and subcritical regimes. In this way, we can 
understand several features of the FF model in terms of 
the same concepts developed for sandpile model. 

The MF equations can be derived by the single site 
approximation to the master equation. Since the deriva- 
tion proceeds as in section IV, we do not repeat it here 
in detail and we present instead some general consider- 
ations. Active sites become stable (fire — > empty) with 
unitary rate, and critical sites become active if ignited by 
the lightning with probability /. The interaction term is 
then given by the fire spreading; an active site create 
as many new active sites as the number of n.n. critical 
sites. To the first order in pa, this term is proportional 
to PaPc times the usual geometrical factor g that takes 
into account the lattice coordination and other model de- 
pendent geometrical effects. The reaction rate equations 
then reads 



Fa 



-pa + fPc + gPcPa + 0{pl). 



(60) 



This expression is very similar to the one obtained in the 
sandpile case, with the exception of the dissipative term 
that here is missing. The two models differ in the dynam- 
ical evolution of stable sites. In the FF model there is a 
term, due of the field p, corresponding to the transition 
rate from stable to critical sites and there is no interac- 
tion between active and stable sites. We can then write 
the stationarity equations for the FF model as 



Pa = fPc + gPcPa 
Pa = PPs 
Pa ^ - Ps ~ Pc, 



(61) 



where we have neglected second order term in pa- As for 
the sandpile model, g is an independent parameter of the 
model and / and p represent the tunable external driving 
fields. The lowest order solutions in / and p to the above 
equations are 



Pa = ^P+^J + 0{P,p') 
Pc =l-li+0{pj) 
Ps =^ + U+0{pJ). 



(62) 



These results have been already obtained in Ref. p7[ , 
where a random neighbor version of the FF is analyzed. 
Their method and the present MF scheme are equiva- 
lent and we will recover the same stationary densities. 
We compute the critical exponents by using the same 
lines adopted for sandpiles, and obtain some new insight 
on the critical properties of the FF model. The den- 
sity of active sites depend linearly upon / and p, which 
are independent driving fields playing the same role as h 
in sandpile automata. If we consider the density of ac- 
tive sites as the order parameter, it appears immediately 
that the critical point is reached if / ^ and p ^ 
simultaneously. This double limit again corresponds to 
the locality breaking of the dynamical rules. In this case 
the order parameter is identically zero in the steady state 
and the system develop long range correlation properties. 
Also for the FF model we can then distinguish among a 
subcritical and a supercritical regime depending on the 
values of the driving fields. 



A. The subcritical regime 

The subcritical and critical regimes correspond to the 
limit in which we have zero order parameter and therefore 
/ = and p = 0. This limit is however not completely 
defined because the density of critical and stable sites de- 
pend upon the ratio f /p. In order to study the critical 
behavior in this limit we repeat the discussion inspired 
by the study of CA with adsorbing states that we already 
used for the sandpiles. Since when / = and p = the 
dynamics is frozen, we have to prepare the system in a 
stationary state in the limit p — )■ and / — > 0, and then 
study the spreading of small perturbations. This is what 
is actually done in numerical simulations where the fire 
evolution and the action of / and p act separately. In 
doing that however, we prepare the system in one of the 
"natural" configurations, corresponding to the stationary 
state in the limit of infinitesimal driving. In this config- 
urations the density of critical sites reaches a limit value 
Pc = 1/5 — 9^^f/p, which depends, via the parameter 
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6 = //p, on the way the limit has been performed. In 
order to study the scahng behavior, we consider the Umit 
/ « p « 9 << 1, keeping 9 constant. In this regime 
we can consider f = p = 0, pa = and we can study the 
system to the leading order in 9. By considering small 
deviations Sp^it) from the stationary state and retaining 
just first order terms in 9, we find the linearized dynam- 
ical equation in diagonal form 



-<5p,(0 = -9Spait). 



(63) 



Hence, the relaxation behavior follows an exponential law 
in which the characteristic relaxation time is given by 
tc ~ 9^'^ with A = 1. This implicitly tells us that the 
system indeed reacts in avalanches. In fact, both driving 
timescales p~^ and are in this regime much longer 
than the characteristic spreading time of an avalanche 
9^^, that therefore remains an isolated event connected 
in space and time. 

Along the lines we followed for sandpiles, we define the 
response function of the system Xf,e{x — x';t — t') that 
characterizes the way the system respond to an external 
perturbation. The response is now function of / and 9 
which is fixed. The total susceptibility Xffi is related to 
the derivative of the stationary density of critical sites, 
and the zero field susceptibility can be obtained as 



Xe 



lim 



Of ■ 



(64) 



Since the density of active sites can be written as 
Pa{f,G) = f/g + {g - l)f/{g9), the singular part of the 
susceptibility diverges as 



xe 



(65) 



In appendix A it is shown that the zero field susceptibil- 
ity is related to the divergence of the average fire size as 
< s >~ X9- Hence, the characteristic fire size is diverging 
for 9 —> 0. This implies that the system is in a subcritical 
regime and perturbations to the stationary state show a 
finite characteristic length for any 9 > 0. Only in the 
limit 9^0 the system responds on all length scales to 
infinitesimal perturbations. We can define the standard 
scaling laws xe = 9^^ with 7 = 1, and ^ ~ 9~'^ that 
characterize the divergence of the correlation length. 

Next, we consider the total response at position r given 
by Xei"!^) = / Xe{''":t)dt. Wc uotc that fire clusters are 
given by the connected clusters of critical sites, because 
in this regime fires are not overlapping. Since a tree can 
burn just once, the average response at distance r is the 
given by the pair connectedness function which is sup- 



posed to behave as 



r(r/0 in MF theory fm. In 



general, by integrating the local response function, we 
have 



(66) 



and therefore by comparing with Eq. ( |65| ) we get 1/ = 1/2. 
It is worth to remark that in this case the above MF rela- 
tions are not enforced by conservation laws and anoma- 
lous exponents can appear in low dimensions. 

To study the avalanche behavior, we introduce the 
probability P{s,9) = s^'^Q{s/ Sc{9)) that a fire involves 
s sites and we identify the usual set of critical exponents 



defined by the scaling laws Sc ~ 9' 



and 



tc ~ . Associated to them we have the scaling relations 
Da — 1/v and 7(7 = 2 — r. We have shown previously 
that in this regime the class of "natural" configurations 
have a density of critical sites which depends on 9, thus 
we consider the difference of densities with respect to the 
critical state 



Pc - Pc{9) - 9^ 



(67) 



We can then find another scaling relation that links the 
avalanche exponents to C, noting that the avalanche size 
distribution corresponds to the distribution of connected 
critical sites cluster. In appendix B we derive this scaling 
relation which results to be 



c = 



T- 1 



(68) 



Collecting all the results obtained above, we have the 
complete set of MF exponents: 



7 = 1, 



3/2, 



D = 4, 



y = 1/2, 



1/2 



(69) 



(70) 



Also in this case, as previously shown by several authors 
|p7| , the MF values correspond to those of mean-field 
percolation. It is important to stress again that in the 
FF case, the absence of a conservation constraint implies 
that MF values for critical exponents are not valid in low 
dimensional system. Anomalous scaling appears below 
the upper critical dimension and the model shows non- 
trivial values of exponents |M . 



B. The supercritical regime 

We consider here the scaling behavior in the region in 
which the order parameter is not zero. In order to re- 
main in the critical region, we must have 9 << 1, but 
now we consider non vanishing / and with / much 
smaller than p. This essentially corresponds to the FF 
model without ignition that in this perspective can be 
considered as the supercritical regime close to the criti- 
cal point. In this limit we obtain immediately from the 
solution of Eq. (|6^) that the order parameter is positive 
and scales as 



Pa ^p'^, 



(71) 



with P — 1. To calculate the relaxation properties we 
have to perform a linear stability analysis of the system 
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around the stationary solutions (|62|), retaining only the 
lowest order terms in p. We consider small fluctuations 
6pK.{t) and the eigenvalues of the diagonal form of the 
dynamical evolutions are 

A± = -{9/2)p±i{gp)^/'. (72) 

The negative real part identifies the characteristic re- 
laxation time that scales as ^ p^^. Together with 
the exponential relaxation, the system shows oscillations 
with period T ^ p^^l"^^ related to the imaginary part 
of the eigenvalues. This MF behavior has been already 
discussed in Ref. 

In the supercritical state the timescale of a perturba- 
tion is comparable to the driving scale, both being of the 
order of Thus, active sites do not spread just on con- 
nected cluster of critical sites. In other words the critical 
sites configuration is not frozen during the perturbation 
and the time evolution connects several clusters of crit- 
ical sites because connecting sites might appear during 
the time evolution. Also in this case, though, the sus- 
ceptibility is given by the total response to a localized 
fluctuation 

Xp~y ^Pait)d±-^V-^- (73) 

Since the response of the system is due to the connectivity 
properties, we have still the usual MF relation Xv ^ ^ ^ 
which implies that ^ ~ p~^ with v' = 1/2. Another way 
to see this result is to think that fluctuation spreads as 
waves of active sites. Since the propagation velocity is 
finite, the correlation length is proportional to the wave 
period T. This simple MF picture does not work in low 
dimensions [Q. 

VI. DISCUSSION AND OPEN QUESTIONS 

A. Relations with Branching processes 

A clear mean-field description of the avalanches in 
SOC models has been obtained through the mapping to 
branching processes p9|-|3^. A branching process is 
defined by a number of active sites that can either die 
or generate n new sites with certain probabilities. The 
simpler example is the case n = 2: a site dies with prob- 
ability 1 — g or generates two new sites with probability 
q. The process usually starts with a single active site 
and continues until no more active sites are present. De- 
pending on the value of q the branching process will die 
after a finite number of steps or continue forever. There 
is a critical value q = qc that separates the two regimes 
(qc = 1/2 for n = 2). For q < qc, the size distribution of 
the branching process is a power law 

P{s) ^ s-3/2/(s/so), (74) 

where the cutoff sq diverges for q = qc- 



It has been shown in Ref. that the Manna model 
can be exactly mapped into a branching process with a 
time dependent parameter q{t) depending on the density 
of critical sites (pc) and on the dissipation A critical 
branching process was obtained as a stationary state in 
the limit of slow driving and conservation |^l| . 

Branching processes can be considered as a general 
framework to describe avalanches in mean-field theory. 
In general terms, we can describe an avalanche by an 
evolving front that can either propagate or stop. In the 
mean-field description, the elements of the fronts do not 
interact and evolve independently. Thus the avalanche 
can be described as a branching process with an effec- 
tive parameter q that depends on the detail of the model 
under study. 

In our formalism, a branching process is associated 
with the propagation of active sites in the subcritical 
regime. In the stationary state for /i = 0, an active site 
generates k = 1, . . . g new active sites with probabilities 

qu^{l-e)(^l^pl{l-pcY-\ (75) 

while no active sites are generated with probability 

go = e + (l-e)(l-Pc)^. (76) 

In this case, the control parameter for the branching pro- 
cess is given by g = kq^ with a critical value = 1- 
In the stationary state, we find pc — 1/g and hence 
q — 1 — e. The critical branching process corresponds 
therefore to the limit e — > 0. A similar analysis can be 
done for the FF model. 



B. Locality Breaking 

We have seen that criticality in stochastic SOC sys- 
tems is achieved only in the limit of infinitesimal driving 
corresponding to the locality breaking of the dynamical 
rules. The non-locality is evident if we consider the zero 
driving limit that is naturally implemented in computer 
simulations using two different timescales, one for the 
avalanche evolution and one for the driving. With this 
infinite timescale separation we introduce new perturba- 
tions only when the system is quiescent; i.e. the evo- 
lution of each site depends on the entire system. For 
a more concrete physical explanation of how the locality 
breaking generates long-range interactions in the system, 
let us consider the case of a vanishing driving rate, corre- 
sponding to a small density of active sites. Because of the 
infinitesimal driving each region devoid of active particle 
is virtually frozen until an active site is generated. The 
activity spreads and in general alter the configuration be- 
fore it moves away or disappears. The active sites leave 
a trace of their dynamical history in the frozen configu- 
rations of critical and stable sites they produced. If new 
active sites are created in the same region at some later 
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times, they can feel the effect of the active sites present 
ear her in the region. This is basically a memory effect, 
which creates a long-range interaction in time and space 
among diffusing active sites. The range of this interaction 
depends on the characteristic timescale of the driving, be- 
cause the fluctuations induced by the driving destroy the 
memory effect. Close to the infinite timescale separation, 
the characteristic driving timescale is diverging and the 
range of the nonlocal interaction extends to the entire 
system. A local interaction is recovered, however, if we 
introduce a size cut-off in the wandering region of active 
particles. This is the case of dissipative sandpilc in which 
after a finite number of steps the active sites disappear 
Q . Over this characteristic size active particles do not 
interact and to obtain a long-range non-local interaction 
the dissipation should go to zero. The same discussion 
applies to the FF model, due to the finite range of con- 
nected critical sites obtained by timing the ratio of / and 
P- 

In this framework, the infinite timescale separation in- 
troduces a nontrivial long-range interactions between ac- 
tive sites which lead to a singular behavior; i.e. criti- 
cal properties in time and space. In sandpile and FF 
model, dissipations or other driving rates introduce a fi- 
nite cut-off to this non-local interactions, which induces 
a subcritical regime. The critical point is reached just 
in correspondence of a second limit corresponding to the 
celebrated double timescale separation. This framework 
leaves room for the appearance of systems in which just 
the single timescale separation could be enough to get 
criticality. These system might show a phase diagram 
without subcritical phase. 



C. Conclusions 

In this paper we have presented a unified mean-field 
theory for stochastic SOC models. We have treated these 
model in analogy with other non-equilibrium cellular au- 
tomata, using a single-site approximation to the master 
equation. With the present approach, we are able to 
identify the order parameter and the control parameters 
of the models and to emphasize similarity and differences 
between SOC and other non-equilibrium system. In par- 
ticular, the language of cellular automata with absorb- 
ing state can be employed to describe SOC models. For 
finite driving rates, we find a supercritical regime char- 
acterized by a finite fraction of active sites. In the limit 
of infinitesimal driving, the system is subcritical and dis- 
plays avalanche response. Criticality arise from a double 
limit: the driving rate and the dissipation (in the sand- 
pile model) or the two driving rates (in the FF model) 
should have vanishing values. This limit correspond to 
the onset of non-local dynamical rules, which are respon- 
sible for the critical behavior characteristic of SOC. 

In this perspective, SOC models appear to be a sub- 
set of non-equilibrium system with steady states, reach- 



ing criticality by the fine tuning of control parameters. 
While this statement is technically correct, we note that 
SOC system are quite peculiar, since the fine tuning can 
only be achieved by limit procedure. This is in contrast 
with ordinary critical phenomena, where the control pa- 
rameter can be directly tuned to its critical value. In 
this sense, SOC systems are less sensitive to fine tuning 
|lO| . Moreover, the driving rate can in general be small in 
many natural phenomena, and this could make the SOC 
framework relevant. 
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APPENDIX A: RESPONSE FUNCTION 
PROPERTIES 

For small perturbations aromid the stationary state, 
the spontaneous microscopic dynamics can be repre- 
sented by introducing the response function. We first 
consider the sandpile case. If we apply a time-dependent 
perturbation h{x, t) to the stationary state, the density 
of active sites changes as 

Apaix,t) = 

XhA^ -x';t~ t')Ah{x',t')d'^x'dt' + OiiAhf) (Al) 

where Xh,e(x — x';t — t') is the response or generalized 
susceptibility function. Here we assume a stationary and 
homogeneous system, i.e. the two point averages depend 
just on the time or space displacement. The above ex- 
pression is valid in the linear regime, only for small varia- 
tions of the perturbing field. We next derive some simple 
properties of the response function for systems whose dy- 
namics is characterized by avalanches. We first consider 
an impulsive disturbance Ah{x',t') = S{t)S'^{x). This is 
a very small perturbation with respect to the total en- 
ergy input J = J h{x)d'^x. In practice it corresponds 
to the addiction of a energy grain on top of the station- 
ary average driving field. Inserting this perturbation into 
Apa{x,t) yields 

Apaix,t)=XhAx-^t)- (A2) 
We ten define the total susceptibility of the system 

Xh,e = f dt [ XhAx;t)d''x, (A3) 
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which quantifies the total response of the system to an 
impulsive disturbance. The total number of active sites 
due to the perturbation is 



(A4) 



Na= dt Apa{x,t)d''x = Xh, 



In absence of external field h ^ 0, the only active sites 
present in the system are due to the delta-perturbation. 
That is, all the active sites are casually connected in space 
and time, thus forming an avalanche whose average size 
is < s >= Na- This is precisely stated by the following 
expression 



Xe 



= < S >, 



(A5) 



which defines a relation between the average avalanche 
size and the zero field susceptibility. As we have seen 
in the previous sections the above expression is at the 
basis of several scaling relations and it explains together 
with conservation the diffusive behavior of the average 
activity. 

Another way to look at a stationary perturbation or 
equivalently to the variation of the stationary averages is 
the following. We consider a different perturbation 



Ah{x',t') = Ah 



for t' < t, 



(A6) 



corresponding to a uniform driving in space and time. 
By changing the variables of integration to t" = t ~ t' 
and x" = X — x' we obtain 



Apa = Ah 



XHA^"\t")d''x"dt" 



(A7) 



V Jo 



Hence, the density fluctuations are time and space inde- 
pendent, as it must be in the new stationary state with 
h ^ h + Ah. Performing the double integral in the right 
term we get the total susceptibility. Therefore 



Apa = AhXh,e, 



(A8) 



from which we obtain that in the stationary state and for 
infinitesimal perturbations 



Xh. 



lim ^ , 

Ah^O Ah 



Apa dpa{h) 



dh 



(A9) 



Notice that from this equation we are able to provide 
a relation between the total response function and the 
divergence of avalanche size 



< s>^Xe 



lim 



dpajh) 
dh ' 



(AlO) 



Eq. AlO states that the zero-field susceptibility in the 



stationary state and the average avalanche size have the 
same singular behavior in the thermodynamic limit. 

We next consider the Forest-Fire model. In this case 
we have the two driving fields / and p and the response 
function depends upon them. The interesting subcritical 



regime is the one in which we take the limit / — > and 
p — > with 9 = f /p << 1. We study the response of the 
system for small perturbation A/ and a fixed value of 6. 
The general expression that characterizes the response of 
the system is given by 



Apa{x,t) = 

XfA^ -x';t~ t')Af{x', t')d^x'dt' + 0((A/)2). (All) 



As for the the sandpile case we can apply a delta- 
perturbation. It follows that, simply rewriting what we 
derived in the sandpile case, we obtain 



Xe 



:< S >, 



(A12) 



where < s > in this case is the average size of fire events. 
In the same way we can consider a stationary pertur- 
bation Af{x',t') — Af for t' < t and by repeating the 
above arguments we recover 



Xf.e 



df ' 



(A13) 



from which follows that the divergence of the average 
size of fires is related to the zero field susceptibility in 
the usual way. 



APPENDIX B: SCALING RELATIONS 

Here we obtain two scaling relations whose derivation 
is straightforward but rather lengthy. 



1. Sandpile model 

Let us consider the fiow decays of activity in the sub- 
critical regime. We define Pa{s,t) the space integrated 
response of an avalanche of s sites. If we assume scaling 
behavior we have that 



Pa{s,t)~s''w(t/t,) 



(Bl) 



where ts is the upper characteristic time of an avalanche 
of s sites and scales as ts ~ s^/^. By imposing the condi- 
tion that J Pq(s, = s we obtain g = l — z/D. We have 
also that w{Q) — 'w{\) ~ and that pa{s,t) is indepen- 
dent of s for small t. This implies that w{x) x^^^^l^ 
as X ^ 0. The total response function is the average of 
the various possible avalanche response 



(B2) 



Pa{t) = / Pais,t)P{s)ds 



that gives after the proper substitution by means of scal- 
ing relations the expression: 



Pait) ~ exp[-(tAc)("-'^/'^"1 



(B3) 



18 



In the MF picture the above relation is consistent with 
the results obtained from the dynamical equations only 
if 



(r-1) 



1 



(B4) 



thus recovering the relations used in Sec. IV. An analo- 
gous results was already obtained in the paper by Tang 
and Bak ||3^. Again we stress that this is not a general 
scaling relations, but an exponents equality valid just in 
MF theory 



Forest-fire model 

The density of critical sites in the stationary configu- 
ration approaches the critical value for 6* as a power 
law 



Ap, = p,{e = 0) - pM ^ 9^. 



(B5) 



In the subcritical regime we have a complete timescale 
separation. Therefore each spreading of activity involves 
just clusters of connected critical sites. This is because 
the tree growth timescale is much longer than the 
activity timescale, thus preventing that new critical sites 
change the connectivity properties of the configuration. 
In this conditions, the probability P(s, 0) to have an 
avalanche of size s scales as the distribution n{s) of con- 
nected clusters with s critical sites times the size of the 
cluster s. This factor takes into account the probabil- 
ity that the ignition process starts on any site of a clus- 
ter of size s. On the other hand, the density of crit- 
ical sites, leaving apart normalization factors, is given 
by Pc ~ / sn{s)ds that is given by the integral of the 
avalanche distribution. We can therefore write 



(B6) 



where we used the explicit form of the avalanche distribu- 
tion. Noticing that g{s/sc{0)) ~ for s > Sc we obtain 
that the main contribution to the above integral is given 

by 



(B7) 



Apc ~ / s '^ds, 
or as a result of the integration 



(B8) 



By substituting Sc ~ 6^^^'^ in the above exp ression and 
requiring the scaling consistency with Eq. ( B5 ) , we finally 
obtain the scaling relation 
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FIG. 2. The density of active site in the BTW model with 
boundary dissipation as a function of the driving rate h is 
plotted for different system sizes L. 



FIG. 1. Phase diagram of a generic sandpile model. We 
include negative values of dissipation corresponding to a net 
addiction of energy during the avalanches. The trivial super- 
critical regime is given by a saturation of the system, which is 
receiving more energy than it can dissipate. The interesting 
region h < e is discussed in the text. 
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FIG. 3. The susceptibility Xx = dpa/dh as a function of 
the dissipation e, for a system with periodic boundary condi- 
tions and size L=64. The line corresponds to the theoretical 
prediction Xe = l/£- 
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